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Abstract 

Reasons for the rising clinical impact of the bacterium Enterococcus faecium include the species' rapid acquisition of adaptive genetic 
elements. Here, we focused on the impact of recombination on the evolution of E. faecium. We used the recently developed 
BratNextGen algorithm to detect recombinant regions in the core genome of 34 £ faecium strains, including three newly sequenced 
clinical strains. Recombination was found to have a significant impact on the E. faecium genome: of the original 1 .2 million positions in 
the core genome, 0.5 million were predicted to have been affected by recombination in at least one strain. Importantly, strains in one 
of the two major £ faecium clades (clade B), which contains most of the £ faecium human gut commensals, formed the most 
important reservoir for donating foreign DNA to the second major £ faecium clade (clade A), which contains most of the clinical 
isolates. Also, several genomic regions were found to mainly recombine in specific hospital-associated £ faecium strains. One of these 
regions (the epa-like locus) likely encodes the biosynthesis of cell wall polysaccharides. These findings suggest a crucial role for 
recombination in the emergence of £ faecium as a successful hospital-associated pathogen. 

Key words: BratNextGen, comparative genomics, phylogenomics, whole-genome sequencing, nosocomial pathogen, 
antibiotic resistance. 



Introduction 

Enterococci are Gram-positive bacteria that commonly inhabit 
the gastrointestinal-tract of healthy humans and other ani- 
mals. However, since the 1980s, the clinical impact of entero- 
cocci has continuously increased and they now rank as one of 
the leading causes of nosocomial infections of the blood- 
stream, urinary tract, surgical wounds, and other sites (Arias 
and Murray 2012; Gilmore et al. 2013). In the past, entero- 
coccal hospital-acquired infections (HAIs) were predominantly 
caused by £ faecalis, but since the 1990s, HAIs have been 
increasingly associated with £ faecium. Nowadays, £ faecium 
HAIs are almost as common as HAIs caused by £ faecalis 
(Hidron et al. 2008). 

An important cause for the rise of £ faecium as a nosoco- 
mial pathogen has been the species' capacity to easily acquire 
novel adaptive traits, including genetic elements encoding 
antibiotic resistance determinants. For example, nosocomial 
infections caused by ampicillin-resistant £ faecium were first 



reported in the United States in the 1980s, after which they 
spread rapidly around the world (Galloway-Peha et al. 2009). 
Currently, ampicillin resistance is reported in more than 80% 
of clinical £ faecium isolates worldwide (Hidron et al. 2008; 
Willems and van Schaik 2009). Other genetic elements that 
have been acquired by £ faecium include esp and other genes 
located on the esp pathogenicity island ICEEfml, the putative 
hyaluronidase gene hyl, and several genes encoding surface 
proteins (Rice et al. 2003; Vankerckhoven et al. 2004; Coque 
et al. 2005; Klare et al. 2005; Camargo et al. 2006; Hendrickx 
et al. 2007; van Schaik et al. 2010), a putative phosphotrans- 
ferase system that contributes to intestinal colonization during 
antibiotic treatment (Zhang et al. 201 3), and specific insertion 
elements, especially IS 7 6 (Leavis et al. 2007). The cumulative 
acquisition of these, but probably many more, genetic ele- 
ments by £ faecium has resulted in the evolutionary develop- 
ment of a genetically distinct £ faecium subpopulation with a 
progressively increasing fitness in hospitalized patients 
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(Willems et al. 201 1). This subpopulation was initially desig- 
nated lineage C1 (Homan et al. 2002) and later renamed to 
clonal complex 1 7 (CC1 7) on the basis of multilocus sequence 
typing (MLST) data analysis (Willems et al. 2005, 2011). 
Phylogenomic studies, using only a limited number of £ fae- 
cium strains, strengthened these previous findings as they 
revealed an ancient phylogenetic split in the £ faecium pop- 
ulation, essentially dividing animal and human clinical isolates 
into one clade (clade A) and human commensal isolates into 
another clade (clade B) (Galloway-Peha et al. 2012; Palmer 
et al. 2012). However, in a recent study using Bayesian- 
based population genetic modeling of an MLST data set of 
more than 1,700 isolates, Willems et al. (2012) found that 
nosocomial £ faecium strains clustered into three, instead of 
one, subgroups, suggesting different evolutionary trajectories 
for modern clinical isolates. 

Pan-genome analyses of £ faecium indicated that £ fae- 
cium has an open genome which means that the total avail- 
able gene pool within this species is essentially unlimited (van 
Schaik et al. 2010). This also suggests that the £ faecium 
genome is highly plastic and few barriers exist for the acqui- 
sition of foreign genetic elements. Evidence for high levels of 
recombination was found in recent comparative genomics 
analyses which noted the existence of hybrid £ faecium 
strains, of which the genomes were an amalgam of clade A 
and clade B backgrounds (Galloway-Peha et al. 2012; Palmer 
et al. 201 2). Similarly, MLST data analysis showed that certain 
£ faecium subpopulations displayed high levels of admixture, 
indicating that up to 14% of the MLST gene sequences orig- 
inated from other £ faecium subpopulations (Willems et al. 
2012). Furthermore, the majority of gene tree topologies of 
individual MLST genes were incongruent (Willems et al. 2005), 
which together with the above suggest that the impact of 
recombination on the population structure and evolution of 
£ faecium may be considerable. However, recombination in 
£ faecium has so far only been studied using either a limited 
set of strains or only few genes (i.e., seven core genes in the 
case of MLST data). Recent sequencing efforts have signifi- 
cantly increased the number of £ faecium whole-genome 
sequences (WGS) available in the public domain. These 
WGS provide an excellent opportunity to study recombination 
on a high-resolution whole-genome level. Here, we use the 
Bayesian Recombination Tracker (BratNextGen) algorithm 
(Marttinen et al. 2008, 2012; Castillo-Ramirez et al. 2012; 
McNally et al. 2013) to study the impact of recombination 
in £ faecium using a set of 34 WGS, including those of 
three clinical strains that were sequenced for the purpose of 
this study. We show that most of the recombination detected 
in clade A originated from a clade B type background and that 
purging recombinant signals from a genome-wide alignment 
prior to phylogenetic tree-building widened the ancient split 
between clade A and B. Moreover, we found that contempo- 
rary clinical isolates share specific recombination signals 



pointing to recombination-driven adaptations to selection 
pressures associated with hospitalization (e.g., antibiotic 
treatment). 

Materials and Methods 

Sequence Information and Genome Sequencing 

WGS were retrieved from GenBank (ftp://ftp.ncbi.nih.gov/ge 
nomes/Bacteria, last accessed August 3, 2013 and ftp://ftp. 
ncbi.nih.gov/genomes/Bacteria_DRAFT, last accessed August 
3, 2013). Enterococcus faecium WGS were downloaded on 
May 2012 and included 31 available genomes: 1 complete 
and 30 draft genomes. Available WGS of 66 other entero- 
cocci were downloaded in July 2012. Besides using publicly 
available WGS, three £ faecium isolates were newly 
sequenced for the purpose of this study. These included 
strains E155, E525, and E1165, which all belong to MLST 
ST17 (Homan et al. 2002; Willems et al. 2005) and have 
been isolated from hospitalized patients in the United States 
(Chicago), Australia (Melbourne), and Italy (Genoa), respec- 
tively (table 1). Genomic DNA was isolated from cell pellets 
of these strains using the Wizard Genomic DNA Purification kit 
(Promega, Leiden, The Netherlands) according to the manu- 
facturer's instructions. The strains were sequenced using 
lllumina GAII sequencing technology generating 50-bp 
paired-end reads from a library with an average insert size 
of 250 bp (Baseclear B.V., Leiden, The Netherlands). After 
confirming the quality of the raw lllumina reads with 
FastQC vO.10.1 (Simon Andrews, Babraham Bioinformatics, 
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/, 
last accessed August 3, 201 3), reads were assembled de novo 
using Velvet v1 .2.01 (Zerbino and Birney 2008). For each 
Velvet assembly, a range of different /c-mer (hash) lengths 
(1 5—49 bp) was empirically tested to gain an optimal N50. In 
addition, Velvet was run with the options "-exp_cov" (ex- 
pected /c-mer coverage) and "-cov_cutoff" (coverage 
cutoff). The expected /c-mer coverage (Ck) was calculated 
using [Ck = Cx (L-k+ 1 )//_], where L and k are the read 
and /c-mer lengths (in bp), respectively. C is the estimated se- 
quence coverage given by [C= L x N/S], where N is the 
number of reads and S is the estimated size of the genome. 
For the three genome assemblies, an estimated £ faecium 
genome size (S) of 2.75 Mb was used. The coverage cutoff 
was set at 0.25 x Ck as suggested by the Velvet developers 
(Zerbino DR, personal communication). 

Functional Annotation and Core Genome Prediction 

The Prokaryotic Annotation pipeline of ISGA (default settings) 
(Hemmerich et al. 2010) was used to annotate the three 
newly sequenced £ faecium genomes, as well as one of the 
available draft genomes (LCT-EF90) for which no annotation 
was available from GenBank. Annotated proteins were 
assigned to Clusters of Orthologous Groups (COG) (Tatusov 
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Table 1 



Whole-Genome Sequenced £ faecium Strains Available in July 2012 



Strain 


Source 3 


Country 


Year 


Reported 
Resistances 6 
and Virulence Genes 


MLSr 


Gen Bank 
BioProject Codes 


E155 


Clinical isolate (faeces) d 


USA 


1995 


AMP, VAN; esp+ 


17 


PRJNA1 92879 


E525 


Clinical isolate (wound) 


AUS 


1998 


AMP, VAN; esp+ 


17 


PRJNA1 92893 


E1165 


Clinical isolate (wound) 


ITA 


1997 


AMP; esp+ 


17 


PRJNA1 92894 


Aus0004 


Clinical isolate (blood) 


AUS 


1998 


VAN; esp-f- 


17 


PRJNA87025 


C68 


Clinical isolate (faeces) d 


USA 


1996 


AMP, VAN; esp+ hyl+ 


16 


PRJNA40855 


TX0082 


Clinical isolate (blood) 


USA 


1999 


AMP, VAN 


17 


PRJNA61137 


E1162 


Clinical isolate (blood) 


FRA 


1997 


AMP; esp+ 


17 


PRJNA47013 


TX0133A 


Clinical isolate (blood) 


USA 


2006 


VAN; esp+ 


17 


PRJNA61139 


TX0133a01 


Clinical isolate (blood) 


USA 


2006 


VAN; esp+ hyl+ 


17 


PRJNA61143 


TX0133a04 


Clinical isolate (blood) 


USA 


2006 


esp+ 


17 


PRJNA61129 


TX0133B 


Clinical isolate (blood) 


USA 


2006 


esp+ 


17 


PRJNA61141 


TX0133C 


Clinical isolate (blood) 


USA 


2006 


VAN; esp+ 


17 


PRJNA61131 


1,231,410 


Clinical isolate (skin and soft tissue) 


USA 


2005 


AMP, VAN; esp+ hyl+ 


17 


PRJNA55719 


DO (TX16) 


Clinical isolate (blood) 


USA 


1992 


AMP; hyl+ 


18 


PRJNA54089 


1,231,502 


Clinical isolate (blood) 


USA 


2005 


AMP, VAN; esp+ hyl+ 


203 


PRJNA55713 


U0317 


Clinical isolate (urine) 


NLD 


2005 


AMP; esp+ hyl+ 


78 


PRJNA47349 


1,230,933 


Clinical isolate (wound) 


USA 


2005 


AMP, VAN; esp+ hyl+ 


18 


PRJNA55701 


1,231,408 


Clinical isolate (blood) 


USA 


2005 


AMP 


582 


PRJNA55721 


E4453 


Commensal isolate (dog faeces) 


NLD 


2008 


AMP 


192 


PRJNA1 79597 


E1071 


Commensal isolate (faeces) d 


NLD 


2000 


VAN 


32 


PRJNA47015 


E4452 


Commensal isolate (dog faeces) 


NLD 


2008 


AMP 


266 


PRJNA179613 


E1679 


Clinical isolate (vascular catheter) 


BRA 


1998 


AMP, VAN; esp+ 


114 


PRJNA47347 


E1636 


Clinical isolate (blood) 


NLD 


1961 


AMP 


106 


PRJNA47345 


D344SRF 


Spontaneous mutant of strain D344R 


USA 


7 


— 


25 


PRJNA46237 


TC 6 


Transconjugant e 


USA 


? 


hyl+ 


25 


PRJNA41101 


E1039 


Commensal isolate (faeces) 


NLD 


1998 




42 


PRJNA4701 1 


1,231,501 


Clinical isolate (blood) 


USA 


2005 




52 


PRJNA55715 


LCT-EF90 


Derived from E. faecium type strain* 


SWI 


? 




76 


PRJNA141665 


E980 


Commensal isolate (faeces) 


NLD 


1998 




94 


PRJNA47017 


Com 15 


Commensal isolate (faeces) 


USA 


2007 




583 


PRJNA55725 


1,141,733 


Clinical isolate (blood) 


USA 


2005 




52 


PRJNA55717 


PC4.1 


Commensal isolate (faeces) 


AUS 


2008 




720 


PRJNA46979 


Com 12 


Commensal isolate (faeces) 


USA 


2006 




107 


PRJNA55723 


TX1330 


Commensal isolate (faeces) 


USA 


1994 




107 


PRJNA55481 



a AII isolates were of human origin unless stated otherwise. 

kwe only list reported ampicillin (AMP) and vancomycin (VAN) resistances here. 

C MLST data were extracted from the MLST website (Imperial College, London; http://efaecium.mlst.net/, last accessed August 3, 2013) or were determined using MLST 
v1.6 (Larsen et al. 2012). 

isolates collected from the faeces of hospitalized patients. Strains E155 and C68 are regarded as clinical isolates because they were representative clones of hospital 
outbreaks. Strain C68 is described in Carias et al. (1998). 

Transconjugant of mating between strains C68 and D344SRF. 

f This strain was derived from the E. faecium type strain and was cultured at 15 and 37°C for more than 4 weeks prior to DNA isolation and sequencing. The E. faecium 
type strain was first deposited to culture collections in 1946, but the impact of this extensive subculturing of this derivative has not been characterized. 



et al. 2001) as described previously (Snel et al. 2002). Protein 
domains were predicted using SMART v7.0 (Letunic et al. 
2012) and Pfam v26.0 (Finn et al. 2010). To determine 
groups of orthologous proteins, the 97,972 annotated pro- 
teins present in the 34 £ faecium strains were used as input 
for an all-versus-all sequence similarity search using BLASTP 
v2.2.24 (Altschul et al. 1997) (default settings, except for: 
-F "m S"; -e 1.0E-05; -z 97,972 [i.e., the total number of 
proteins]). From the BLAST output, groups of orthologous 



proteins were predicted using OrthoMCL v2.0.2 (Li et al. 
2003). Orthologous groups with exactly one representative 
protein from each of the 34 input strains were considered 
to be part of the £ faecium core genome. 

Detection of Recombination 

The BratNextGen algorithm (Marttinen et al. 2008, 2012) was 
used to detect recombination events in £ faecium. A 
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concatenated multiple sequence alignment of the £ faecium 
core genome was used as input for this algorithm. This con- 
catenated alignment was built as follows: for each ortholo- 
gous group, the corresponding nucleotide sequences were 
extracted and aligned using Muscle v3.7 (Edgar 2004), after 
which gaps were stripped from the alignment using trimAI 
v1.2 (Capella-Gutierrez et al. 2009). Stripped alignments 
were concatenated in the order of orthologous group num- 
bers as assigned by OrthoMCL. Because OrthoMCL assigns 
orthologous group numbers on a "first-come, first-served" 
basis, orthologous genes ended up in the concatenated align- 
ment in the order in which they were read by OrthoMCL, thus 
following the order in which they are located on the genomic 
contigs and scaffolds. As the first genomic contigs that 
entered OrthoMCL were those of strain 1,141,733, the align- 
ment followed the exact gene order in the contigs of this 
strain. Analysis of the alignment showed that the gene 
order used also exactly matched that of 1 1 other strains and 
nearly exactly matched that of the remaining 22 strains (a 
median of only 3 synteny breaks was observed for these 22 
strains with a maximum of only 10 and 1 1 synteny breaks for 
strains Aus0004 and 1,231,408, respectively). Thus, the con- 
catenated alignment of orthologous genes accurately reflects 
the synteny of the average E. faecium genome. The estimation 
of recombination was performed with the default 
BratNextGen settings as in (Marttinen et al. 2012) starting 
out with a partition of the genome alignment into 5kb 
blocks and using 20 iterations of the estimation algorithm, 
which was assessed to be sufficient since changes in the 
hidden Markov model parameters were already negligible 
over the last 30% of the iterations. The alpha parameter 
was estimated using the default options in BratNextGen and 
the resulting value was 4.96. Significance of a recombining 
region was determined as in (Marttinen et al. 2012) using a 
permutation test with 100 permutations executed in parallel 
on a cluster computer (threshold of 5% was used to conclude 
significance for each region). 

Phylogeny and the Identification of Recombinant 
Sequence Origins 

Phylogenetic trees were built using RAxML v7.2.8 (Stamatakis 
2006) under the GTRCAT model. To assess the effect of 
recombination on the evolution of E. faecium, two phyloge- 
netic trees were built: one using the aligned SNPs from the 
concatenated core genome alignment, and one using the 
same information, but filtered for significant recombination 
signals as determined by BratNextGen. Confidence was in- 
ferred by running 1,000 bootstrap replicates under the 
GTRCAT model. The phylogenetic origin of the identified 
recombinant sequences was assessed as follows: first, for 
each recombination event, the corresponding sequence align- 
ment was extracted from the larger E. faecium core genome 
alignment and a separate phylogenetic subtree was built. 



Recombinant sequences were classified as belonging to the 
same recombination event when they had received exactly 
the same recombination start and end positions by 
BratNextGen. To prevent obscuring the individual subtrees 
with signals from other recombination events, all sequences 
that partially or completely belonged to another recombination 
event were removed from the alignment prior to building the 
tree (note that recombination events do not overlap within one 
sequence). Second, the resulting trees were midpoint rooted 
using the Phangorn package v1.6-5 (http://cran.r-project.org/ 
web/packages/phangorn/index.html, last accessed August 3, 
2013) in R (R Development Core Team, http://www.R-proj- 
ect.org, last accessed August 3, 2013), after which a decision 
scheme was followed using the topology of the first major split 
closest to the root (supplementary fig. S1, Supplementary 
Material online). In the whole-genome recombination-filtered 
tree, E. faecium was separated into the two clades A and B (see 
also Results and Discussion). This knowledge was taken into 
account in the decision scheme, such that recombinant se- 
quences were predicted to have originated from either an (an- 
cestral) clade A or clade B strain. Alternatively, when the 
recombinant sequences clustered separately from both 
clades (i.e., as an outgroup), the sequence was predicted to 
have originated from a species other than E. faecium. Because 
BratNextGen uses a conservative approach for predicting 
recombinant sequences to avoid false positives (Marttinen 
et al. 2012), it is likely that not all recombination events were 
detected. Therefore, undetected overlapping recombination 
events could still obscure individual subtrees and blur the sep- 
aration into clades A and B. To take this into account, a max- 
imum of two supposedly nonrecombinant sequences were 
deleted from the tree if needed to get a "pure" clade A/B 
separation, thereby keeping a minimum of three clade A and 
three clade B strains in the tree. If the tree could not be purified 
in this way, the phylogenetic source of the recombinant se- 
quence was considered to be either E. faecium or inconclusive 
(supplementary fig. S1, Supplementary Material online). 

Analysis of Recombinant Sequences Originating from 
Outside £ faecium 

In cases where recombinant sequences were predicted to orig- 
inate from outside £ faecium, an attempt was made to iden- 
tify the potential donor species. For this, recombinant 
sequences were split up into their corresponding gene 
sequences (partial gene sequences <30 bp were not consid- 
ered). Recombinant gene sequences were then screened 
against GenBank's nucleotide sequence database (nt) (ftp:// 
ftp.ncbi.nlm.nih.gov/blast/db/, last accessed August 3, 2013; 
date of download August 201 2) and against the assemblies of 
66 non-£ faecium enterococci (discussed earlier) using 
BLASTN v2.2.24 (default settings, except for: -F F). Target se- 
quences were recovered using an F-value cutoff of 1.0E-10 
and a query sequence coverage of >75%. Recovered 
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sequences were aligned with the corresponding E. faecium 
gene sequences, after which midpoint rooted individual 
gene trees were built using RAxML. Trees were purified as 
described above if necessary and were inspected manually 
to assign potential donor species. 

Database Submission 

The whole genome sequence data of E. faecium strains E1 55, 
E525, and E1165 have been deposited at DDBJ/EMBL/ 
GenBank under the accessions AUWXOOOOOOOO, 
AUWV00000000, and AUWW00000000, respectively. 

Results and Discussion 

Enterococcus faecium Phylogeny Based on the Core 
Genome 

To study the effects of recombination on the genome evolu- 
tion of E. faecium, we used the genome sequence information 
of 34 E. faecium strains, including 31 publicly available strains 
and 3 strains that were newly sequenced for the purpose of 
this study. The newly sequenced strains were E1 55, E525, and 
E1 165 and were all modern clinical isolates, originating from 
hospitalized patients in the United States, Australia, and Italy, 
respectively, and belonging to ST17, which is an important 
node in the hospital-associated subpopulation of E. faecium 
(Homan et al. 2002; Willems et al. 2005, 2011, 2012) 
(table 1). De novo assembly of the sequence reads resulted 
in an average genome size of 2.94Mb (±0.1 4 Mb) and a 
mean scaffold N50 of 43.4 kb (±2.9 kb). The number of scaf- 
folds of size >200 bp ranged from 1 85 to 242 (supplementary 
table S1, Supplementary Material online). Annotation of the 
new genomes using the ISGA pipeline resulted in an average 
set of 2,971 (±170) protein-encoding genes per genome. 

A total of 97,927 protein sequences were extracted from 
the 34 E. faecium genomes and used as input for assessing the 
core genome. This resulted in a set of 5,075 orthologous 
groups, of which 1,512 were considered to belong to the 
core genome as they contained exactly one protein represen- 
tative from each of the 34 input strains. After stripping the 
individual core gene alignments for gaps, all alignments were 
concatenated to generate a whole-genome core alignment of 
1,223,137 bp. From this alignment, we extracted 85,488 
SNPs, which were used for the construction of a phylogenetic 
tree. The topology of the resulting tree (fig. 1/4 and supple- 
mentary fig. S2, Supplementary Material online, for details on 
bootstrap support) was in agreement with what has been 
observed before by others who used an E. faecium core 
genome extracted from 8 to 29 isolates to build a phyloge- 
netic tree (Galloway-Peha etal. 2012; Lam etal. 2012; Palmer 
et al. 2012; Willems et al. 2012). The tree displayed a deep 
phylogenetic split between two major E. faecium clades. 
These clades were previously observed by Palmer et al. 
(2012), who designated them clades A and B, and by 



Galloway-Peha et al. (2012), who designated them clades 
HA (hospital-associated) and CA (commensal/community-as- 
sociated), respectively. Clade A contained most of the clinical 
and hospital-associated isolates as well as animal isolates, 
whereas clade B consisted mainly of human commensal iso- 
lates. The clinical isolate 1,231,408 clustered consistently 
(100% bootstrap support) with clade A, albeit relatively dis- 
tantly from all other clade A strains. This is also in line with 
previously built phylogenetic trees and with the observation 
that 1 ,231 ,408 is a hybrid strain that appears to have received 
a large part of its genome from a clade B background (Palmer 
et al. 201 2) (discussed later). As expected on the basis of their 
ST (table 1), the newly sequenced clinical isolates clustered 
closely together with other modern clinical isolates of the 
same ST, as well as with strain C68 (ST16), which is a single 
locus variant of ST17 (Homan et al. 2002). 

Enterococcus faecium Phylogeny Based on a 
Recombination-Filtered Core Genome 

To estimate the effects of recombination on the evolution of 
E. faecium, we used the algorithm BratNextGen (Marttinen 
et al. 2008, 2012). Running this algorithm on the E. faecium 
core genome alignment resulted in the prediction of 1 .2 Mb of 
recombinant sequence, divided over 383 recombination 
events (supplementary fig. S3, Supplementary Material 
online). Of the original 1,223,137 positions in the core 
genome alignment, 534,961 (44%) were predicted to be 
affected by recombination in at least one strain. On the SNP 
level, this rate was comparable with 38,665 of the original 
85,488 positions (45%) identified as being recombinant. 
The remaining 46,823 nonrecombinant SNPs were used to 
build a recombination-free phylogenetic tree. As shown in 
figure 1B (supplementary fig. S2, Supplementary Material 
online, for details on bootstrap support), the topology of 
this tree was largely similar to that of the unfiltered tree pre- 
sented in figure 1/4 in the sense that the ancient split between 
clade A and B was still observed after purging for recombina- 
tion. Interestingly, BratNextGen did not detect any recombi- 
nation signals in clade B. This does not necessarily mean that 
clade B strains are completely unaffected by recombination. 
BratNextGen uses a conservative approach for predicting re- 
combinant sequences to avoid false positives (Marttinen et al. 
2012). Consequently, true recombination signals may go 
undetected (i.e., false negatives) when there is low statistical 
power. The reasons for the absence of detectable recombina- 
tion in clade B may be that 1) clade B strains were underrep- 
resented in the available WGS data set used here, and 2) the 
diversity within clade B was relatively large as compared to the 
diversity within clade A (fig. 14). Because BratNextGen has 
higher statistical power to detect significant recombination 
events when they are shared across strains, these factors 
(i.e., limited sample of strains and high diversity between 
them) may have hampered the detection of recombinant 
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Fig. 1. — Impact of recombination on £ faecium phylogeny and quantification and determination of origins of recombinant sequences. (A) Midpoint 
rooted phylogenetic tree (RAxML) built from an alignment of 85,488 variable positions in the E. faecium core genome. No recombination filtering was 
applied. (6) Midpoint rooted phylogenetic tree (RAxML) build from an alignment of 46,823 nonrecombinant variable positions in the E. faecium core 
genome. Numbers near branches indicate percentage bootstrap support: only numbers <95% are indicated and only for main branches (for all bootstrap 
supports, see supplementary fig. S2, Supplementary Material online). Branch lengths correspond to the scale bar, in units of changes/nucleotide position. 
(0 Predicted levels and origin of recombinant sequences per strain. 



sequence segments within clade B. Indeed, we noticed that, 
when using our decision scheme for predicting recombinant 
sequence origins (supplementary fig. S1, Supplementary 
Material online), tree purifications were made 74 times corre- 
sponding to the removal of 103 sequences. This observation 
may point to a recombinant nature for the corresponding se- 
quences although other reasons for sequence divergence (i.e., 
mutation) cannot be ruled out. Of the 103 removed 
sequences, 60 belonged to clade B strains, mostly to strains 
LCT-EF90 (33 cases) and E980 (23 cases), thus potentially in- 
dicating the presence of recombinant sequences in clade 
B. Nevertheless, the fact that BratNextGen detected no re- 
combination signals at all in clade B at least strongly suggests 
that, in general, clade A is more prone to recombination than 
clade B. Consequently, the recombination-filtering procedure 
effectively meant a purification of clade A only. As became 
evident by visually comparing the trees in figure 1/4 and B, the 
purification resulted in a strongly increased cohesiveness of 
the phylogenetic tree for clade A. To quantitatively support 
this observation, we extracted all strain-to-strain pairwise 
branch-length distances from both trees and compared the 
average distance within clade A and B and between both 



clades, pre- and post-recombination-filtering. Potential over- 
sampling of lineages was taken into account by including only 
one strain (the one furthest from the root) for branches where 
all pairwise strain distances were less than 5.0E-04 in the 
purified tree (e.g., strain TX0133a01 was picked as a repre- 
sentative for all 5 TX0133 strains). This analysis showed that 
upon recombination-filtering, the average distance within 
clade A significantly decreased (from 0.1 to 0.03; P<0.01; 
Wilcoxon signed rank test), whereas the average distance 
between clade A and B significantly increased (from 1 .4 to 
2.0; P<0.01). The average distance within clade B (0.2) 
remained the same. In summary, recombination-filtering wid- 
ened the ancient clade A/B split by bringing the clade A strains 
closer together. 

Comparison of the two trees showed that most of the 
strains that clustered tightly together in the unfiltered tree 
also clustered tightly together upon recombination-filtering. 
The most notable exceptions were 1) the separation of strains 
E1039 and 1,231,501 from each other and from the larger 
cluster containing for instance strain E1 679, 2) the shift of the 
5 x TX0133 branch to a position within the ST17 cluster of 
modern clinical isolates, and 3) the shift of the branch 
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1,231,502/110317 to a position just outside the ST17 cluster 
(fig. 1 and supplementary fig. S2 [Supplementary Material 
online] for details). These last two modifications resulted in a 
robust branch consisting only of modern clinical isolates of 
ST17 (and 16) and with ST78 and ST203 isolates together 
forming a distinct lineage. Another major change that 
occurred upon recombination-filtering was the displacement 
of strain 1,231,408 from a "hybrid" position in-between 
clade A and B toward a position branching deeper inside 
clade A and tightly clustering with strain 1,230,933, which 
clustered with strain DO in the unfiltered tree. The strong 
shift of strain 1,231,408 toward clade A and away from 
clade B fits with the previously observed hybrid nature of 
this strain (Palmer et al. 201 2). However, a similar post-recom- 
bination-filtering displacement of strain 1,231,408 was not 
observed in a previous study (Qin et al. 2012), in which the 
authors used PHI (Bruen et al. 2006) to detect and remove 
recombination signals prior to phylogenetic tree building. In 
their resulting tree, strain 1,231,408 still occupied a "hybrid" 
position similar to the one observed in the unfiltered tree 
presented here (fig. 1/4). 

Quantities of Recombination in Clade A 

Recombination signals were detected in each of the 27 clade 
A strains included in this study. However, the amount of de- 
tected recombinant DNA as a percentage of the 1 .2 Mb core 
genome varied widely over the different clade A strains 
(fig. 10- When correcting for potential over-sampling of ge- 
netic lineages, the average amount of recombination detected 
in clade A was 4.4% of the core genome. The strain with the 
highest amount of recombination was strain 1,231,408 with 
as much as 26.9% of its core genome predicted to be recom- 
binant. Other strains with a relatively high load of recombina- 
tion were E1 039 (1 0. 1 %) and 1 ,23 1 ,501 (8.2%). At the other 
end of the spectrum were strains TX0133a01, E1162, and 
TX0082 with only 0.3%, 0.7%, and 0.7% of their core 
genome predicted to be of foreign origin, respectively 
(fig. 10- Interestingly, the modern clinical isolates of ST16 
and ST17 had a significantly lower amount of recombination 
(average of 1.1%) than the other strains in clade A (6.0%) 
(P < 0.01 ; Mann-Whitney U test). This finding appears to be in 
agreement with a previous observation made by Willems et al. 
(2012), who applied a genetic admixture analysis on MLST 
sequence data from a large collection of £ faecium strains 
and found that hospital isolates displayed lower levels of re- 
combination than isolates from nonhospitalized persons and 
pigs. This may point to genetic and/or ecological isolation of £ 
faecium hospital isolates and fits with the previously proposed 
hypothesis that novel, highly invasive, £ faecium strains may 
arise through horizontal gene transfer, but once they have 
adopted their new lifestyle, they become genetically isolated 
from other bacteria and become less prone to recombination 
(Willems et al. 2012). 



The Recombinant DNA in Clade A Mainly Had a 
Clade B Origin 

Interestingly, the vast majority of the £ faecium recombinant 
sequence was predicted to have originated from a clade B 
type £ faecium background (fig. 10- More specifically, we 
found that of the 1 .2 Mb of recombinant sequence, 71 % had 
a predicted clade B type source, followed by 19.6% with a 
predicted source from outside £ faecium and 7.6% with a 
predicted clade A type origin. For only 0.1 % of the recombi- 
nant sequence, it was unclear whether it originated from a 
clade A or clade B type background. Also, for only 1 .7% of the 
recombinant sequence, the phylogenetic source was incon- 
clusive according to the decision scheme that we employed 
(supplementary fig. S1, Supplementary Material online). These 
results indicated that clade B £ faecium strains have formed, 
and possibly still form, an important source for donating core 
genomic DNA to clade A £ faecium strains. Especially strains 
1,231,408, E1039, and 1,231,501 were found to have 
received a large amount of core genomic DNA from a clade 
B type source, totalling 327, 103 and 84 kb, respectively. Also 
with respect to the total amount of foreign DNA detected in 
their core genomes, the amount of DNA with a potential clade 
B type source in these strains was high at 99%, 84%, and 
84%, respectively. In figure 2, we plotted the detected recom- 
binant sequences and their potential sources onto the core 
genome sequence of the completed strain Aus0004, as 
determined in this study. The results presented in figure 2 
expand on previous results of Palmer et al. (2012), who stud- 
ied genome mosaicism in five clade A and three clade B £ 
faecium strains. For each gene in these eight strains they de- 
termined whether it had either a clade A- or clade B-type 
signature. Mosaicism was most evident for two clade A strains 
they analyzed, namely strains 1,231,408 and 1,231,501. 
Here, when using the currently much larger data set of 34 
£ faecium genomes, we found that strains 1,231,408 and 
1,231,501 still ranked among the top strains with respect to 
the amount of recombinant sequence coming from a clade B 
type background. In accordance with the observations made 
by Palmer et al. (2012), we found that these recombinant 
sequences clustered into a few long genomic stretches 
taking up around 27% and 7.8% of the 1,231,408 and 
1,231,501 core genomes, respectively (fig. 2). 

Recombinant DNA with a Predicted Source Outside 
£ faecium 

As mentioned earlier, almost 20% of the recombinant core 
genomic DNA in £ faecium was predicted to come from a 
source outside £ faecium. To find potential donor species for 
these recombinant segments, we extracted the corresponding 
gene sequences from each recombinant segment and 
screened them against a database of enterococcal genome 
sequences (Materials and Methods) and GenBank's nt data- 
base. Similar sequences were extracted and used, together 
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with the corresponding recombinant and nonrecombinant E. 
faecium gene sequences, to build phylogenetic trees 
(Materials and Methods). The resulting trees were manually 
inspected to find close homologs of (i.e., potential donor spe- 
cies for) the recombinant E. faecium sequences. In cases 
where one recombination event covered multiple genes, all 
the resulting gene trees had to be in agreement with each 
other concerning potential donor species. If not, for instance 
because one of the recombinant genes did not yield BLAST 
hits outside E. faecium, the donor species was deemed un- 
known. Unfortunately, we did not find potential donor species 
for the majority (84 out of 1 1 1) of the corresponding recom- 
bination events. This was due to the fact that for these 84 
events, no similar sequences were detected for at least one of 
the recombinant genes, other than E. faecium sequences 
themselves. The most likely explanation for obtaining so few 
BLAST hits may be that current microbial sequence databases 
simply do not (yet) contain sequences from the actual donor 
species or from species closely related to these donors. Still, for 
the remaining 27 recombination events (24%), we did find 
more specific evidence concerning the nature of the potential 
donor species. For instance, for 25 out of the 27 recombina- 
tion events, the recombinant gene sequences clustered in be- 
tween £ faecium and other enterococcal species or clustered 



together with multiple different enterococcal species other 
than E. faecium, suggesting that these recombinant se- 
quences have originated from the genus Enterococcus (results 
not shown). Finally, for two recombination events, we found 
even more specific evidence about the nature of the potential 
donor species. One of these events was found in strain C68 
and involved genes efxg_00252-efxg_00259. Each of these 
genes clustered together with a gene from E. hirae ATCC 
9790 (ehr_06405-ehr_06440, respectively), suggesting that 
E. hirae or a closely related species has been the responsible 
recombination donor for these genes. Similarly, we found that 
gene efme1071_1769 in strain E1071 had an E. faecalis sig- 
nature. These observations indicate that there is also detect- 
able recombination between E. faecium and other species in 
the genus Enterococcus thereby significantly enlarging the 
potential reservoir for genetic diversity in E. faecium. 

Functional Analysis of Recombinant Genes 

To describe the functional impact of recombination on the 
E. faecium genome, we analyzed whether certain functional 
gene categories were more affected by recombination than 
others. To this aim, we assigned COGs to each orthologous 
protein group in the E. faecium genome. We then compared 
the composition of the COG functional categories found in 
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the recombinant regions versus those found in the 
recombination-free regions. After correcting for potential 
over-sampling, we found that 1,832 proteins were affected 
by recombination compared with 40,502 recombination-free 
proteins, corresponding to 2,086 and 45,654 assigned COG 
functional categories, respectively. As shown in supplemen- 
tary figure S44, Supplementary Material online, a linear 
regression analysis between the COG category compositions 
of these two data sets showed a strong linear correlation 
(/? 2 = 0.87), suggesting that in general most functional 
groups are equally prone to recombination in the £ faecium 
core genome. Still, when we used Grubb's test for outliers 
using a 99% confidence level, we found that COG categories 
C (energy production and conversion) and S (function un- 
known) were identified as outliers, potentially being relatively 
more (category C) and less (category S) prone to recombina- 
tion in the E. faecium core genome (supplementary fig. S44, 
Supplementary Material online). Further analysis of the cate- 
gories C and S showed that their skewed representation in 
recombinant DNA was not due to specific COGs that were 
over- or underrepresented, but was caused by a proportional 
change of all COGs in the entire category (results not shown). 
By a similar approach, we also evaluated whether specific 
functional categories originated from specific recombination 
donors. Supplementary figure S46-D, Supplementary 
Material online, shows the COG distributions for the recom- 
binant genes originating from one of the three major donors 
of recombinant DNA (clade A, clade B, and outside £ faecium) 
compared with the COG distributions for the recombinant 
genes originating from all donors, except the one of interest. 
This analysis revealed that recombinant genes belonging to 
COG category G (carbohydrate transport and metabolism) 
mainly originated from clade A. Two individual COGs from 
category G that were overrepresented in recombinant DNA 
from clade A, were COG2723 and COG1486, which both 
encode putative beta-glucosidases. Together these two 
COGs made up 62% of all recombinant COGs of category 
G originating from clade A, whereas they only made up 9% of 
all recombinant COGs of this category originating from other 
sources. Other functions that originated mostly from clade A 
included a putative sugar uptake permease (COG4975), a 
dihydroxyacetone kinase (COG2376), and an ABC type poly- 
saccharide transporter (COG 1 134) (results not shown). These 
results suggest that genes involved in sugar uptake and 
metabolism are extensively transferred between clade A 
£ faecium strains. This is of interest because genes involved 
in these processes may be important for successful coloniza- 
tion of the gut during antibiotic treatment, as has been dem- 
onstrated recently (Zhang et al. 2013). 

Functional Analysis of Recombination Hotpsots 

To locate the position of conserved major recombinant 
regions, we searched for recombination events that were 



detected in at least two strains (correcting for potential 
over-sampling) with a total summed length of at least 5 kb. 
This revealed a set of 18 recombination events (listed in sup- 
plementary table S2, Supplementary Material online), which 
merged into eight regions of overlapping and closely flanking 
recombination events (fig. 2). Interestingly, six of these eight 
regions (regions 2-7), including the four largest ones with a 
length ranging from 26 to 59 kb, were located in the third 
quadrant of the Aus0004 core genome. These data indicate a 
nonrandom distribution of recombination events over the 
£ faecium core genome and suggest the existence of a hot- 
spot for recombination in £ faecium, which, when translating 
to the Aus0004 genome, roughly lies between open reading 
frames efau004_01 354 and efau004_02178. Besides identi- 
fying recombination events that were conserved across mul- 
tiple isolates, we searched for sites in the £ faecium core 
genome with clear traces of recombination-driven gene ac- 
quisition. For this, we screened for genomic regions of at least 
five consecutive accessory genes that were flanked on both 
sides by the same recombination event. The accessory genes 
and the flanking recombinant core genes had to be located on 
one single contig or scaffold. A total of 17 recombination 
events met these criteria, several of which were clearly asso- 
ciated with mobile genetic elements, such as phages, trans- 
posons, recombinases, and integrases (supplementary table 
S2, Supplementary Material online). Of the 17 events, five 
were found in strain 1,231,408, corresponding to the first 
quadrant of the Aus0004 core genome. Still, 9 of the remain- 
ing 12 events were found within the recombination hotspot 
described earlier, which further indicates that this region on 
the £ faecium genome is relatively prone to recombination. 

Recombination of an Ampicillin Resistance Determinant 

As mentioned earlier, we found that the clinical £ faecium 
strains of ST1 6/1 7 displayed relatively low levels of recombi- 
nation. Considering this, it is of interest to note that 6 of the 
18 conserved recombination events listed earlier were associ- 
ated with ST16/17 strains, indicating that the recombination- 
driven diversity of these genomic regions may reflect 
adaptation to selection pressures associated with hospitaliza- 
tion. Four of the corresponding recombination events (events 
8-11, supplementary table S2, Supplementary Material 
online) were located within region 5 (fig. 2). Among these 
were events 8 and 9, which were predicted by BratNextGen 
to be independent recombination events, but which covered 
exactly the same set of core genes. Together, these two 
events were found in all ST16/17 strains analyzed, except 
TX0133. Other events that partly overlapped with events 8 
and 9 were found in strains E4452 and E1071 (results not 
shown). One of the genes within the affected region was 
ddcP (efau004_01934 in Aus0004), which has previously 
been identified as a membrane-associated D-alanyl-D-alanine 
carboxypeptidase that confers resistance to ampicillin and 
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lysozyme in £ faecium E1162 (Zhang et al. 2012). 
Recombinant ddcP was predicted to stem from an ancestral 
clade B type £ faecium. On the protein level, three positions 
gave recombinant DdcP a clear clade B type signature, includ- 
ing an lle-to-Val at position 1 01 , an Asp-to-Ala at position 298 
and an Asp-to-Glu at position 392 (supplementary fig. S5, 
Supplementary Material online). Although these amino acid 
variations do not by themselves fully explain DdcP-mediated 
ampicillin resistance in £ faecium (e.g., strain E1071 and all 
clade B strains are ampicillin sensitive), they do suggest a sce- 
nario in which clinical strains, like E1162, have recombined 
their ddcP gene with a clade B type ddcP gene, which in 
combination with other factors like mutations in pbp5, has 
resulted in clinical strains becoming more prone to developing 
ampicillin resistance. 

Recombination of the epa-Like Locus Is Reminiscent of 
Pneumococcal "Capsule-Switching" 

The two remaining recombination events that were associated 
with ST16/17 £ faecium strains were both found in region 2 
(fig. 2) and were predicted to have originated from an uni- 
dentified source outside £ faecium. Both events covered and 
flanked almost the same set of core and accessory genes, 
respectively, namely efau004_01354 to efau004J385 in 
strain Aus0004 and its orthologs in strains E1 165 and E525 
(event 2, supplementary table S2, Supplementary Material 
online) and a somewhat larger region extending toward the 
orthologs of efau004_01387 in strains 1,231,410, E155, 
1,230,933, DO, 1,231,502, and U0317 (event 3). In 
Aus0004, the region between efau004_01354 to 
efau004_1387 mostly consists of accessory genes, including 
two ISL3-family transposases and a retroviral integrase, further 
pointing toward a history of lateral gene transfer for this com- 
plete genomic region. The region involved has previously been 
described as a downstream extension of the £ faecium epa- 
like locus (Palmer et al. 2012; Qin et al. 2012), a gene cluster 
that is putatively involved in the biosynthesis of antigenic cell 
wall polysaccharide, based on homology with the £ faecalis 
epa locus (Xu et al. 1998; Teng et al. 2009). The extended 
region of the epa-like locus (corresponding to efau004_ 
01365-efau004_0l382, plus efau004_01 384 in Aus0004) is 
known to be highly variable across £ faecium and has been 
divided into four locus variants (Qin et al. 2012). Analysis of 
the current set of 34 £ faecium strains showed that the same 
division into four locus variants still holds. Figure 3 shows rep- 
resentatives of these four locus variants as well as the regions 
affected by recombination. Notably, all the epa-like loci af- 
fected by recombination events 2 and 3 belonged to locus 
variant four. Only one epa-like locus of this type, the one 
present in the clade B strain E980, did not contain recombi- 
nation signals according to BratNextGen. This does not nec- 
essarily mean that this locus has not recombined in this strain, 
but may suggest that the potential recombination event was 



not recent enough, so that any strong recombination signals 
may have been neutralized over time. Besides recombination 
events 2 and 3, additional independent events were detected 
within the same genomic region of strains E1165, E4453, 
E1679, E4452, E1636, and E1039, but these events covered 
only small segments within epaMNO, never extending beyond 
three consecutive core genes (results not shown). 
Nevertheless, the observation that multiple independent 
recombination events were detected in the downstream 
epa-like extension, suggests that this region of the £ faecium 
genome is prone to recombination. As shown in figure 3, the 
variable epa-like region that was flanked and thus potentially 
affected by events 2 and 3 contained a predicted class C beta- 
lactamase gene (efau004_01 384 in Aus0004) and sialic acid 
biosynthesis genes neuABCD (efau004_01 366 and 
efau004_01370-efau004_01372) specific for locus variant 
four. Sialic acid decoration of the bacterial cell surface may 
be a form of molecular mimicry that affects detection by the 
host's immune system (Vimr and Lichtensteiger 2002). 
Interestingly, events 2 and 3 both extended beyond the epa- 
like locus into a core genomic region containing yet another 
class C beta-lactamase gene (efau004_01360). A recombi- 
nant cell surface biosynthesis cluster flanked by two beta-lac- 
tamase genes has also been described in Streptococcus 
pneumoniae. In this species, serotype variants have arisen 
through recombinational replacements within and around 
this so-called cps locus (Trzcihski et al. 2004; Brueggemann 
et al. 2007; Croucher et al. 2011; Golubchik et al. 2012; 
Marttinen et al. 2012). Our observations suggest that recom- 
bination may have a role in capsule variation in £ faecium. 

Conclusions 

In this study, we have analyzed the impact of recombination on 
the whole-genome level in £ faecium. We show that filtering 
for recombination prior to phylogenetic tree building had a 
purifying effect on the ancient £ faecium clade A/B split. 
Closer inspection of the origin of the recombinant sequences 
indicated that this purifying effect was mainly caused by the 
fact that clade A strains have received foreign DNA mainly from 
clade B type strains. However, no recombination signals were 
detected in clade B strains themselves, suggesting that the 
transfer of DNA within the £ faecium population is mainly a 
unidirectional event: from clade B to clade A, and not vice versa. 
However, it is important to note that the number of clade B 
strains in our data set was relatively low, which may have ham- 
pered the detection of recombination signals within this clade. 
The inclusion of more clade B strains in future recombination 
detection analyses will identify whether the apparent lack of 
recombination within clade B is genuine or is merely a result of a 
lack of resolution. Still, the results presented here strongly indi- 
cate that clade A is more prone to recombination than clade B. 
Given the observation that recombination seems to act as a 
coalescent force for clades A and B, it is of interest to note 
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Fig. 3. — Enterococcus faecium epa-like locus variants and recombination signals detected in locus variant four. Gene clusters of the four epa-like locus 
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downstream epa-like extension is indicated by blue solid lines above the illustration of locus variant four. Locus variant four was predicted to be affected by 
large core genomic recombination events 2 and 3 (indicated by solid black lines above the locus variant four illustration) that potentially spanned multiple 
accessory genes (indicated by dashed gray lines). Genes are color-coded according to functional category. Drawings are to scale. 



that, within clade A, contemporary clinical isolates of ST16/17 
contained relatively low levels of recombination. This fits with 
previous hypotheses (Willems et al. 2012; McNally et al. 2013) 
that novel, highly successful, clinical strains may arise through 
gene acquisition, but become genetically isolated and less 
prone to recombination once they have adopted their new 
life-style. Still, several recombination events were associated 
with contemporary clinical isolates. These events involved the 
ampicillin resistance determinant ddcP and the putative cell wall 
polysaccharide biosynthesis epa-like locus. These findings point 
to important recombination-driven mechanisms that facilitate 
adaptation of E. faecium to selection pressures associated with 
hospitalization. 

Supplementary Material 

Supplementary figures S1-S5 and tables S1 and S2 are avail- 
able at Genome Biology and Evolution online (http://www. 
gbe.oxfordjournals.org). 
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